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Abstract 


This report presents multigrid methods for solving the 3-D incompressible viscous rotating 
flows in a NASA low-speed centrifugal compressor and a marine propeller 4119. Numerical for- 
mulations are given in both the rotating reference frame and the absolute frame. Comparisons are 
made for the accuracy, efficiency, and robustness between the steady-state scheme and the time-ac- 
curate scheme for simulating viscous rotating flows for complex internal and external flow applica- 
tions. Prospects for further increase in efficiency and accuracy of unsteady time-accurate computa- 
tions are discussed. 



I. Introduction 


Computation of unsteady rotating flows can be applied in many practical areas such as in turbo- 
machinery cascade flow, helicopter rotor flow, and marine propulsor flow predictions[l][2][3]. Two 
numerical approaches are typically used in the simulation of unsteady flows in rotating machinery. 
A common way is to solve the governing equations in a rotating reference frame by computing either 
the relative velocity [4] or the absolute velocity [2][5]. The advantage of this approach is that it uses 
a steady-state formulation, if the flow field can be viewed as a steady state in the rotating frame. 
Thus, many efficient acceleration techniques, such as local time stepping and multigrid method, can 
be used. An alternative to the steady-state approach is to establish the governing equations in a fixed 
absolute frame, and solve the equations using a time-accurate formulation. However, the demand 
on accuracy and efficiency for time-accurate solutions is much higher than that for steady-state 
solutions. For the reason of accuracy, time step is restricted, and it has to be chosen smaller than 
the smallest characteristic scale length to be resolved. The restriction on the time step reduces the 
efficiency of implicit schemes, but it is the most straightforward way to deal with the general un- 
steady flows that can not be viewed as a steady state in the rotating reference frame, such as the un- 
steady rotor/stator interaction. 

This report describes and compares both numerical formulations of solving the governing equa- 
tions in the rotating reference frame with a steady-state method and in the absolute frame with an 
unsteady time-accurate method. In the rotating frame, a steady relative velocity flow field is pur- 
sued. In order to fully use the original code which is written in the absolute velocity components, 
a dependent variable transformation is first performed to change the relative velocity components 
into the absolute velocity components in the governing equations, and the computation is performed 
based on absolute velocity instead of relative velocity in the rotating frame [2] [5]. Two practical 
applications are presented for solving viscous rotating flows in a NASA low-speed centrifugal com- 
pressor and in a marine propeller 4119. The purpose here is to evaluate and validate the accuracy. 
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efficiency and robustness of the current method to predict complex flow fields for both internal and 
external flow applications. 

In the following, the 3-D incompressible Navier-Stokes equations in general curvilinear coordi- 
nates are first given, followed by the numerical methods used to solve the equations in both the rotat- 
ing reference frame and absolute frame. The difference between the two approaches is addressed. 
Then the multigrid implementation [6] [7] is illustrated to accelerate the solutions in both the rotating 
and absolute frames. Computational results of viscous rotating flows in both impeller and propeller 
cases are presented. In the last section, some conclusions are given. 
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n. Governing Equations 


The artificial compressibility form of the 3-D incompressible Navier-Stokes equations in a gen- 
eral curvilinear coordinate system (§, rj, £, x), which rotates about the x-axis at a constant speed of 
Q, can be written as follows 


d Q , dF , dG , dH _ c 

# + f + + 


( 1 ) 


where Q, the flux vectors F,G,H and the source term S are 


Q 



p ' 


0 ' 


u 


0 

= J 

V 

5 = 7 

-Qw 


. w . 


. Qv . 


F = J 


P(U' - £,) 
uU' + % x p - 
VU' + ZyP - f f ' 

wU 1 + ZzP ~ 

■ z 


G = J 


W ~ V t) 
uV' + rj x p - 
vV + rj y p - 
wV + VzP ~ F 

’ z 


H = J 


fKW - tt) 

uW + txp - \ 

VW + lyP - 

wW + t, zP - % 

^2 


3 



Here J, p, u, v, and w denote the Jacobian transformation, pressure, and the Cartesian velocity com- 
ponents in the absolute frame, respectively, is the artificial compressibility coefficient, and terms 
where fc-f , tj, and £, are the viscous flux components in curvilinear coordinates. The 
Baldwin-Lomax algebraic turbulence model is adopted in this work. The relative contravariant ve- 
locity components If, V, and W are defined as 

U' = i x u! + % y v' + w' + 

V' = rj x u' + rj y v' +t] z w'+ rj t 
U' *£,«' + v' + ? z W +Zt 

where u', v 1 , and w' are the relative velocity components in the rotating reference frame, and are writ- 
ten in terms of absolute velocity components u, v, and w as 

u' = u v' = v — Qz w' = w + Qy 

To solve the governing equations in an absolute frame, one can simply set the rotating speed 
Q to zero in the above equations and evaluate £ f , rj t , and & as a result of grid rotation, which reduces 
to the normal conservative form of 3-D incompressible Navier-Stokes equations in general curvili- 
near coordinates based on the fixed absolute frame. Grid speeds have been included in the above 
formulation to allow grid motion relative to both the rotating and the absolute frame. In this work, 
however, only steady state solutions with stationary grids are pursued in the rotating frame. These 
steady state solutions of the absolute velocity on stationary grids in the rotating frame are carried 
out by setting Q to the rotational speed and £/ = ty = & = 0 in Eq. (1). These absolute velocity solu- 
tions , which are viewed as steady state in the rotating frame, correspond to a particular time and posi- 
tion of the rotating grid in the absolute frame. The unsteady computations with dynamic grids are 
performed in the absolute frame. 
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HI. Numerical Solution Method 


In both the steady-state and unsteady formulations, the governing equations (1) are discretized 
by a cell centered finite-volume scheme. The time derivative is differenced using the Euler back- 
ward formula. For the one-dimensional case in both the rotating and the absolute frames (in the 
absolute frame, with the source term S equal to zero), it may be written as 


Q1 


n + 1 


Q n t 


Ax 



= s ? +1 


( 2 ) 


where the index i corresponds to a cell center and indices i ± 1/2 correspond to cell faces. The spatial 
discretization of the Euler fluxes at cell faces are approximated by using the third-order MUSCL 
approach in the Roe scheme. Details about this method can be found in [8]. 

The nonlinear system of equations (2) is solved by the discretized Newton-relaxation (DNR) 
method described in [9]. Note that Eq. (2) can be written in a simple form as 


N (Q n + l ) = S n+l 

Applying Newton’s method to Eq. (3) yields 

N '(Qn + l,m^Qn+l,m+l _ Qn+l,m ^ _ _ (N(Q n+l ’ m ) _ gn+l,m^ 


( 3 ) 


( 4 ) 


where m= 1, 2, 3,... is the number of Newton iterations implemented at each time step, with <2 rt+1 > 1 
=Q n . These sub-iterations at each time level serve to eliminate linearization error, and thereby 
insure temporal accuracy. N 1 is the Jacobian matrix of the nonlinear equation (3) where the con- 
tribution of the source term is not included. The resulting formulation of Eq. (2) by Newton’s 
method is 
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where F 1 is the Jacobian matrix of the numerical flux vector F, with the first subscript representing 
the position of the cell face of the numerical flux vector, and the second subscript representing the 
position of the dependent variable vector that the numerical flux vector is differentiated with re- 
spect to. I a is an identity matrix except the first diagonal element is zero in order to satisfy the true 
incompressible continuity equation. In the rotating frame where the flow field is a steady state, 
continued iteration of Eq. (5) would presumably lead to Q n+l — ► Q n . For unsteady computation in 
the absolute frame, time-dependent results can be achieved when a converged solution is obtained 
at each time step through Newton iterations. Gauss-Seidel relaxations are used to solve the linear 
system of equations, which results from Newton’s method, approximately at each Newton itera- 
tion. 
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IV. Multigrid Method 


In this work, the multigrid method [6] [7] is used to accelerate both the steady-state and time-ac- 
curate computations. The difference between the steady and unsteady multigrid methods is that in 
the former, time is advanced in the fine grid as well as the coarse grid to achieve full efficiency. 
While in the later, both the fine grid and coarse grid equations must be solved at the same time level 
to ensure temporal consistency [10]. The two-level multigrid method for Eq. (3) can be briefly de- 
scribed as follows: 

1 . Iterate N^Q^S 11 X times on the fine grid h by Newton’s method. 

2. Restrict the residual and solution to the coarser grid 2h, and iterate N 2h (Q 2k ) = S 2h J'T 
times, where S 2h =N 2h (I 2 \Q k ) + R 2 \(S h - is the source term on the coarse grid 2h. 

3. Interpolate the correction from the coarser grid to the fine grid and update the solution 

Q h <= Q h +P h 2h(Q 2h ~I 2h hQ k ) 

4. Repeat steps 1~3 for Jib times at the same time level, using Q h as the new approximation 
to Q n+1 . 

In the above procedure, X is the number of Newton iterations for the fine grid and coarser grids, 
and jWj is the number of multigrid cycles implemented at each time step. Choosing different values 
of J'T and Jit "may form different multigrid strategies and result in different effects. In this work, 
the number of Newton iterations JT is chosen as one for all computations. The number of multigrid 
cycles used at each time step varies according to different cases. Note that if no coarser grids 
exist, both parameters Jib and JT have the same meanings in regard to the implementation of the 
code. 
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V. Results 


Two applications involving both internal and external viscous flows are presented here. One is 
a turbulent flow in a NASA low-speed centrifugal compressor (LSCC) [ 1 ] . Another case is a turbu- 
lent flow about a marine propeller 41 19 [3]. In each case, computations were first performed by a 
steady numerical approach in the rotating frame, then compared with the unsteady time-accurate 
approach in the absolute frame. Computations were carried out on a single processor of an SGI 
75MHz R8000 workstation. 

A NASA Low-Speed Centrifugal Compressor 

The study of the LSCC is sponsored by the NASA Lewis Research Center to evaluate the capa- 
bility of the computational method to predict the flow field in the complex geometric channels of 
centrifugal compressors. The complex phenomena to be considered include secondary flows in the 
impeller passage, tip clearance flows, etc. The experimental investigation of the LSCC was con- 
ducted by Hathaway et al. in [1], which gives detail measurements of the velocity components V r> 
Vt (Vq in Ref.[l]), and V x in the blade passage. These experimental data are used to verify the accura- 
cy of computed results. 

Two computational meshes for the LSCC geometry were built with about 300K and 200K grid 
points on the fine and coarse grids, respectively. The grid spacing on the fine grid surfaces is 4.X10 -6 , 
and has 73 points on the blade, 13 points on the tip clearance, 37 points spanwise, and 4 1 points pitch- 
wise (Figure 1). The Reynolds number is 4.3xl0 6 , based on the velocity of the inlet flow and the 
diameter of the blade tip. The grid y + value on surfaces is about 1. The inflow boundary condition 
is given by specifying the three velocity components u, v, and w. The back pressure is specified at 
the exit of the computational domain. Figure 2 shows the convergence history of the steady-state 
solution in the rotating frame, based on the fine grid. A 3-level full coarsening multigrid strategy 
was employed to accelerate the convergence. A local time stepping was used with a maximum CFL 
number of 20. At each time step, only one multigrid cycle (Jti=l), with 5 Gauss-Seidel relaxations 
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was implemented. The residual was reduced by 2 orders of magnitude in 200 multigrid cycles, and 
then became flat. The reason for that may be due to the extremely small volume in this fine grid 
(minimum volume is less than l.xl0~ 14 ). The computer requirement for this solution is 336 MB 
internal memory, and 27 hours of CPU time for 500 multigrid cycles on the machine mentioned 
above. 

The following Figures 3(a) through 3(d) show computed and measured relative velocity compo- 
nents at 4 locations (m/m s = 0.149, 0.475, 0.644, and 0.941, where m/m s is the non-dimensional 
shroud meridional distance) along the blade passage. The results in each plot are shown at every 
5% of span, with the results nearest the shroud located at 95% of blade span from the hub, where 
100% span denotes the blade tip. Agreement between the numerical prediction and the measurement 
is considered reasonably good. Figures 4(a) through (d) show computed particle traces on the suc- 
tion and pressure surfaces, and the top and hub surfaces of the LSCC channels, based on the relative 
velocity. The plots indicate a continuing migration of fluid outward toward the tip near the blade 
surface (Figures 4 (a) and (b)). 

The LSCC flow field was also solved based on the coarse mesh, using both the steady approach 
in the rotating frame and the unsteady time-accurate approach in the absolute frame. For the steady 
computation, the same multigrid strategy and parameters were employed as in the previous case. 
Figure 5 shows the convergence history in the coarse grid. A convergent solution was obtained in 
500 multigrid cycles by reducing the residual for 4 orders of magnitude. The comparison of relative 
velocity components between the computation and experiment is very similar to the fine grid case 
(shown in Figures 3(a)~(d)), and therefore are not presented here. The unsteady calculation was 
started using local time stepping, and moving the entire grid at a rate of 300 time steps per revolution 
of the impeller. A 3-level multigrid procedure was also used to accelerate the convergence of the 
unsteady solution at each time step. After the solution marched for two revolutions, the calculation 
was switched to minimum time stepping, with 3000 time steps per revolution (A £=0.00029). The 
final solution was obtained after about 4000 time iterations, due to the very small time step used in 
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the unsteady computation. The accuracy of the time-dependent solution will be discussed below. 
Another effort made by Taylor [11] using the Roe flux formulation (non-MUSCL type [12]) and 
a single grid approach showed similar results after about 6000 time step iterations of the unsteady 
computation. The cost of CPU time in these unsteady solutions is about 6~ 1 0 times that of the steady 
solution in the rotating frame. 

The accuracy of time-dependent solutions is assessed with different time steps and multigrid 
cycles at each time step, by comparing with the results obtained using a steady approach in the rotat- 
ing frame. The unsteady calculation can be performed at a larger time step of 300 cycles per revolu- 
tion (A 1=0.0029), but the results differ from those obtained using the steady method, see Figure 6(a). 
Applying more multigrid iterations (Jib) at each time step seemed not to improve the solution very 
much, as shown in Figure 6(b). The accuracy of the unsteady solution is significantly improved only 
after reducing the time step to 0.00029, which is 3000 time steps per revolution (Figure 6). The 
above numerical results reveal the following characters of the code. First, the current unsteady code 
is robust, since the allowed time step is not bounded by numerics (stability), but by physics (resolu- 
tion). Therefore, a largest possible time step can be selected in the computation to achieve best effi- 
ciency while maintaining desired accuracy. Second, one multigrid cycle is sufficient to provide a 
convergent solution at each time step in this case, since more iterations do not change the solution 
significantly. Third, the accuracy of time-accurate solutions is critically dependent on the time step 
used in the calculation, which suggests that higher order time accuracy may be preferred in unsteady 
computations, especially in complex internal flows. Though using a very large time step (or local 
time stepping) does not provide sufficient accuracy for the time-dependent solution, it is an efficient 
way to obtain an initial approximation to start with, or to quickly predict the flow field qualitatively. 
This strategy is also adopted in the computation of the next case. 

A Marine Propeller 4119 

There is reported computational result for the marine propeller 4119 [13]. The purpose here is 
to further assess the accuracy and efficiency of the current method for simulating external rotating 
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flows in a complex geometry. The computational grid consists of three blocks with a total of 280K 
points (Fig.7). There are 41 points in the streamwise and spanwise direction on the blade. The Re- 
ynolds number is 5.76x10 s , based on the freestream velocity and the diameter of the blade. Figure 
8 shows the convergence history of the steady solution obtained in the rotating frame, where the cal- 
culation was performed by using a 4—level multigrid full coarsening. At each time step, one multi- 
grid iteration (Jb= 1), with 5 Gauss-Seidel relaxations was implemented. Local time stepping was 
used at a CFL number of 5. It is seen that the residual is reduced by about 3 orders of magnitude 
in 400 multigrid cycles. The CPU time of this calculation is about 20.6 hours on the SGI R8000 
workstation. 

The pressure coefficient distributions on the blade surface are show in Figures 9(a)~(c), where 
r is the radial distance from the measured point on the blade to the hull axis, and R is the radius of 
blade tip. Favorable agreement was obtained between the computation and experiment, except for 
the pressure side at location r/R=0.3. 

A computation was also performed using the unsteady time-accurate approach in the absolute 
frame. The initial solution was obtained by running the unsteady code in the absolute frame using 
local time stepping, while the computational grid was rotated at a rate of 200 cycles per revolution 
of the propulsor. After 200 time steps, the time-accurate calculation was started using a minimum 
time step of 0.004165, which is equivalent to 200 time steps per revolution. The multigrid cycles 
are employed to insure the convergence of the solution at each time step. It was found that the final 
solution with one multigrid cycle (Jlb=I) at each time step is close to that with two multigrid cycles 
(J1 d= 2) at each time step. The unsteady solution became periodic after seven revolutions of the grid 
motion, or 1400 time steps. The cost of the CPU time of this unsteady solution is 3.5 times and 7 
times that of the previous steady solution, with one and two multigrid cycles at each time step, re- 
spectively. Figures 10(a)~(c) show computed w-velocity contours obtained in the rotating frame and 
in the absolute frame with different multigrid cycles. Again, results obtained in the two different 
frames are similar. Figures 1 1 (a)~ 1 1 (c) show the computed and measured pressure coefficient dis- 
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tributions on the blade surface, where computed results were obtained in the absolute frame using 
2 multigrid iterations at each time step. Computational results about thrust and torque coefficients 
are given in the Table 1. A desired accuracy is achieved in both computations. 
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VI. Conclusions 


The computational results are presented for both internal and external viscous turbulent flows 
in complex geometries, a NASA low-speed centrifugal compressor and a marine propeller 4119. 
Computations are performed and compared in both the rotating reference frame using a steady-state 
formulation and the absolute frame using a time-accurate approach. A multigrid strategy is applied 
in both the steady-state and time-accurate computations to improve the efficiency and robustness 
of the algorithm. Results show that solving the unsteady rotating flow in a rotating frame costs much 
less CPU time than solving the flow in the absolute frame, in both internal and external flow cases. 
Comparison between the computed results and the experimental data is considered satisfactory for 
the current turbulence model and grid resolution. 

Although the current multigrid algorithms are relatively efficient and robust, there are some un- 
resolved issues which require further investigation. These include the use of a larger CFL number 
in the steady flow computations with a source term, and selection of the optimal time step for the 
best output regarding both accuracy and efficiency in time-dependent solutions. Finally, the second 
or higher order time difference scheme should be incorporated into the code to improve the accuracy 
and efficiency for time-accurate computations. 
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Figure 3 (a) Comparison of computed (in rotating frame) and measured relative 
velocity magnitude at station 85 (m/m s =0.149) 
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Figure 3 (b) Comparison of computed (in rotating frame) and measured relative 
velocity magnitude at station 118 (m/m s =0.475) 
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Figure 3 (c) Comparison of computed (in rotating frame) and measured relative 
velocity magnitude at station 135 (m/m s =0.644) 
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Figure 3 (d) Comparison of computed (in rotating frame) and measured relative 
velocity magnitude at station 165 (m/m s =0.941) 
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Figure 6 Comparison of computed velocity components at station 118 at 
95% span at different multigrid cycles (JHo) and time steps 
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Figure 10 Computed u-velocity contours (a) in rotating frame, 

(b) in absolute frame with 1 multigrid cycle, 

(c) in absolute frame with 2 multigrid cycles 
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Figure 1 1 Computed (in absolute frame) and measured pressure coefficients 
on the blade surfaces at (a) r/R=0.3, (b) r/R=0.7, (c) r/R=0.9 
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Thrust Coefficient 

Torque Coefficient 

Experimental Data 

0.146 

0.028 

Computation 
in Rotating Frame 

0.1497 

0.0254 

Computation 
in Absolute Frame 

0.1498 

0.0256 


Table 1, Measured and computed thrust coefficient and torque coefficient 















